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SUMMARY 


ADAM is an acronym for Ax i symmetric Duct Aeroacoustic Modeling and is the name 
of an interconnected system of computer programs for analysing the propagation 
and attenuation of sound in aeroengine ducts containing realistic compressible 
subsonic mean flows. It has been developed primarily for research directed 
towards the reduction of noise emitted from turbofan aircraft engines and 
includes as two basic components, a streamtube curvature program for determin- 
ation of the mean flow and a finite element code for solution of the acoustic 
propagation problem. 

The system, which has been specifically tailored for ease of use, is presently 
installed at NASA Langley Research Center on a Control Data Cyber 175 Computer 
under the NOS Operating System employing a Tektronix terminal for interactive 
graphics. 

This report describes the scope and organization of the ADAM System and in- 
cludes a users guide, examples of input data and results for selected cases. 
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1 INTRODUCTION 


BACKGROUND 


ADAM, which is an acronym for Ax i symmetric Duct Aeroacoustic Modeling, has its 
origins in research work on computational duct acoustics sponsored by NASA- 
Langley Research Center since the mid-1970's. This work had as its goal the 
achievement of a better understanding of the physics of sound propagation and 
attenuation in lined aeroengine ducts containing near sonic compressible mean 
flows. 

In performance of this work, tools have been developed as research vehicles 
both for ad hoc use in particular projects, and also with the impetus to ad- 
vance computational algorithms to their limits. 

From this perspective, the ADAM system should be viewed as a byproduct of 
research in aeroengine duct acoustics rather than as the result of a pre- 
determined research or development objective. In this context, the system 
described in this report is an existing capability whose maintainance and 
development is dependent upon Langley's future orientation and goals in duct 
acoustics. 


STATEMENT OF THE PROBLEM 


Mathematically, the linearized, inviscid, duct acoustics problem is usually 
formulated as a boundary value problem in quantities representing small pertur- 
bations of the basic conservation of momentum, energy and continuity equations 
of fluid dynamics. Two distinct sets of equations result from the perturbation 
process and these are considered to represent two distinct phenomena within 
the physical system: 

■ the unperturbed mean aerodynamic flow field, which is itself 
unaffected by the presence of acoustic energy propagating through it 
and 

■ the acoustic field whose total energy content is unaffected by the 
mean flow but whose local energy propagation characteristics are 
affected (sometimes dramatically) by it. 

This assumption is only one of many simplifications and generalizations inherent 
in the process of creating a general purpose computer code to model a class of 
physical problems. Such assumptions, which lie not only in the realm of physics 
but also in mathematics and computer science, are sometimes theoretically justi- 
fiable, but ultimately rely on careful comparisons with laboratory work for 
their verification. 

In the absence of flow, combined research in computational and experimental duct 
acoustics has already been instrumental in verifying the theoretical formulation 
of acoustic boundary conditions * and variable geometry effects 3 . 
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Now, with the evolution of a model such as ADAM, basic comparisons between 
experiment and theory regarding the effects of flow on sound propagation and 
attenuation in ducts can be performed. 


PREVIOUS WORK 


Computational Acoustics, that is the discipline of obtaining numerical solu- 
tions to the partial differential equations of acoustics, had its beginning 
in the 1970‘s and is expected to come of age in the 1980's. 

Linear duct acoustics has been a leading research area in the development of 
this methodology. Baumeister 4,5 presented successive state-of-the-art revues 
on numerical techniques in linear duct acoustics citing 86 references in his 
November 1981 paper 5 . This single statistic indicates the substantial com- 
pendium of literature which has accumulated in this field. 

Four distinct techniques from the work conducted to date may currently be 
viewed as contenders for the title of "Most Favored Method in Linear Duct 
Acoustics" : 

■ Steady state finite element 

■ Steady state finite difference 

■ Transient finite difference, and 

■ Wave envelope method. 

Steady state finite element methods, so called because they rely on a spectral 
decomposition in time, are presently the best developed of these techniques. 

The finite element method offers the advantage of easily coping with complex 
geometries and when combined with a direct matrix solution algorithm yields 
reliable robust computer codes. The principal disadvantage of the method is 
the substantial computational cost for large problems. 

Steady state finite difference methods suffer from the disadvantage of re- 
quiring involved geometry mapping while still consuming computational resources 
comparable to steady state finite element methods. For simple geometries, how- 
ever, the implementation in computer code is simpler than for steady state 
finite elements. 

Transient finite difference techniques have been developed by Baumeister 6,7,8 
with the motivation of reducing computational requirements relative to steady 
state methods. A brief note of explanation regarding the importance of compu- 
tational requirements in evaluating different techniques is necessary. All 
methods referred to above (including wave envelope techniques) require sub- 
stantial computer resources for realistic problems. Because of this fact, the 
magnitude of the differences in computer cost between techniques is highly sig- 
nificant, either for reduced cost, or added capability at comparable cost. 
Iterative solution of sparse matrices where possible is usually faster than 
direct solution. Steady state methods are not amendable to iterative solution 
algorithms because they usually give rise to nonpositive definite or ill-con- 
ditioned matrices 9 . 
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The transient finite difference method, however, is a special case which is 
amenable to iterative solution and holds the potential for future advances in 
the state-of-the-art. It does have the disadvantage however of requiring the 
same involved geometry mapping as the steady state finite difference method and 
requires further research before it can be regarded as robust, even in simple 
problems. 

Wave envelope methods also base their appeal on reduced computing cost but 
through a very different precept. Although several different wave envelope im- 
plementations exist, they all attempt to utilize the physical reality that while 
the basic acoustic variables vary rapidly in space in wavelike patterns, their 
amplitudes do not. Thus, formulations based on amplitudes may be performed on 
a much coarser grid than formulations based on the variables themselves. 

Examples of wave envelope formulations are contained in references 10 through 13. 
Generally these methods continue to hold considerable appeal for solution of 
acoustic problems. 


SYNOPSIS 


Principally because of its generality, advanced state of development, and over- 
all robustness, the acoustic analysis implemented in ADAM is a steady state 
finite element method. This algorithm was originally developed with iso- 
parametric linear elements 14 - 16 and has been modified to include the optional 
use of subparametric Hermite elements. 

The flow model used in ADAM is the Streamtube Curvature Program 17 developed by 
the General Electric Company. It is a well tested and thoroughly documented 
program designed specifically for numerical computations of compressible flow 
about and within aeroengine nacelles. Its operation is based on the Streamtube 
Curvature Relaxation Technique 18 and employs the concept of summing a number of 
psuedo one-dimensional confluent streamtubes to obtain the total flow. To 
append viscous boundary layer effects onto the inviscid model, an iteration is 
performed to match a separate boundary layer calculation over each surface with 
the overall streamtube solution 19 . 

User interface with both codes has been simplified and the connection between 
streamtube curvature flow program and acoustic finite element program has been 
completely automated. Interactive graphics for aerodynamic and acoustic results 
complete the system. 

Since ADAM has evolved from research activity in a young active research field, 
the techniques implemented within it, while considered the best currently avail- 
able, are expected to become dated with time. Nevertheless, it is anticipated 
that the existence of ADAM will be by itself a stimulus to further work and an 
aid to the verification of new techniques. 

Section 2 of this document contains an overview of the system and describes its 
general organization. It also forms an introduction to the users guide con- 
tained in Section 6. Sections 3, 4, and 5 describe the basic components of the 
system namely, the program to determine the mean flow, the acoustic finite 
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element code, and the programs for interactive graphics. Sections 7 and 8 
contain examples and some comments on the current status of verification and 
testing. 
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2 SYSTEM OVERVIEW 


The ADAM system consists of five FORTRAN computer programs interconnected 
through a defined set of data files. It is presently installed at NASA- 
Langley Research Center on a Control Data Cyber 175 computer under the NOS 
operating system. User interface with ADAM has been strongly oriented towards 
Tektronix interactive terminals. 

All FORTRAN code has been specifically designed to be reasonably hardware inde- 
pendent. Major hardware dependence is through the "Tektronix Plot 10" graphics 
system, but this is not considered a major barrier to portability since several 
other terminals emulate the Tektronix. 

The description of the system contained in this report, however, is limited to 
its installation at NASA-Langley and portability is not guaranteed. 

With the aid of the simplified functional diagram contained in Figure 2.1, con- 
sider a step-by-step walkthrough of a typical modeling sequence. 

Phase A: Determination of Mean Flow 


(1) Preparation of the input file for the General Electric Streamtube 
Curvature Program. 

Principal items to be specified are: 

o Free stream Mach number, ambient pressure, temperature, and 
fluid properties. 

o Parameters to control automatic refinement of spatial com- 
putational grid. 

o Duct geometry. 

(2) Execution of GE-STC Program. 

The STC program computes the subsonic and transonic field of 
inviscid compressible flow about and/or within the given geometry. 
An optional boundary layer analysis evaluates the friction losses 
and displacement of the inviscid flow. Since the computational 
mesh does not extend into the boundary layer, its surfaces are dis- 
placed by the computed boundary layer thickness. A one-seventh 
power law Mach number profile is assumed within the boundary layer. 

Within limits prescribed by the user, the program computes its own 
computational grid. Ambient pressure, temperature, and fluid prop- 
erties, as well as duct geometry and program output in the form of 
mean flow parameters at nodes of the computational mesh, and are 
written to an output data file. 
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Phase A: Determination of Mean Flow 



Figure 2.1 Simplified Functional Diagram of an 
ADAM Modeling Sequence 


The user may view streamlines, Mach number contours, or velocity 
contours plotted over localized sections or over the entire analysis 
area, and may obtain hard copies of these plots for a permanent 
record. 

Phase B: Solution of the Acoustic Propagation Problem 

(1) User preparation of the input file for Acoustic Finite Element 
Program. 
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The acoustic finite element program has been arranged to give the 
user maximum flexibility for problem specification while at the same 
time minimizing the quantity of input data required. This has been 
implemented by allowing the user to specify input in the form of 
FORTRAN functions. 

Principal input parameters are as follows: 

■ Number of finite elements in r and z directions 

■ Location of duct section for acoustic analysis within aerodyna- 
mic coordinate reference frame 

■ Frequency and spinning mode number 

■ Type of finite element to be used in the analysis (linear or 
Hermite) 

Input functions which may be constants, data tables, or analytical 
functions, include the following: 

■ Axial and radial grid point distributions 

■ Inner and outer wall acoustic admittance 

■ Termination admittance 

■ Source-plane pressure distribution. 

(2) Execution of AFE Program. 

The AFE program obtains duct geometry from the GE-STC program output 
file and creates the mesh for the acoustic finite element analysis. 

It then obtains the parameters defining the mean aerodynamic flow on 
the GE-STC mesh and interpolates within this mesh to define the flow 
field on the acoustic mesh. 

Utilizing Hermite or linear finite elements as defined by input data, 
the program assembles the global finite element matrix by numerically 
integrating over each finite element and applying boundary conditions 
as specified in user input admittance and source functions. 

A block tridiagonal solver is then applied to obtain primitive 
acoustic variables consisting of pressure and three components of 
velocity. Acoustic intensity is integrated over a radial plane at 
each axial station to give net acoustic flux. All the above results 
are written to a program output data file. 

(3) User Inspection of Results via Interactive Graphics. 

Due to the large number of complex nodal variables and the corres- 
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pondingly greater number of possible interpolation points, some 
form of graphics is essential for proper user inspection of results. 

The interactive graphics program for interrogating the output data 
has been specially developed for this purpose with the capability 
not only of displaying stored output data, but also of interpolating 
within Hermite (or linear) elements or computing special functions 
such as admittance along a boundary. 

The user may request plots of amplitude or phase of any analysis 
variable along any axial or radial nodal line, specifying optional 
plot scales. Hard copies of these plots may also be obtained for a 
permanent record. 



3 THE COMPRESSIBLE FLOW PROGRAM 


GENERAL 


The Streamtube Curvature Analysis (STC) method and its implementation as a com- 
puter program for the analysis of aeroengine nacelle flow field properties and 
drag, is fully described in reference 17. This original program has been 
slightly modified to produce an output file for use by the Acoustic Finite Ele- 
ment (AFE) program. 

A brief summary of the STC theory and computer program structure, together with 
user instructions are included in this report for completeness. For details 
not included here, the reader is referred to the original publication. 

The STC method solves the inviscid equations of motion over a two-dimensional 
body (plane or axisymmetric) at transonic speeds. The boundary layer procedure 
of Stratford and Beavers (SAB) (reference 19) is incorporated to permit evalua- 
tion of boundary layer displacement effects and friction drag. 

The program may be used to compute both internal flow through the engine nacelle 
in a number of flow channels as well as external flow around the nacelle. This 
capability is important for producing accurate flow and drag computations. The 
acoustic finite element program is currently limited, however, to a single axi- 
symmetric or planar section so that a two-stage analysis as described in Section 
6 is necessary for real engines. 

Although the program has the capacity to define and iteratively improve its com- 
putational mesh, the user has control over the amount of refinement and total 
mesh size. Geometry is specified in the form of a coordinate list with optional 
derivatives. Boundary layer effects may be selectively included for any surface 
and imposed at any level of grid refinement. 

THE STREAM TUBE CURVATURE METHOD 

The STC program is designed to solve the equations of motion along streamlines 
[v = constant lines) and along lines which are orthogonal to the streamlines 
(c = constant lines). The variable s is introduced to avoid confusion with the 
velocity potential <t> which is only applicable when the flow is irrotational . 

Across the streamlines, the continuity and Crocco form of the momentum equation 
are written: 


Continuity: 

3A = ” 

ptfl 


(c 

= Const. ) 

(3.1) 

Momentum 

i = -d 2 

5 an v 

9H _ 8S 
an an 

(c 

= Const.) 

(3.2) 
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Along the streamlines the following forms of the energy and momentum equations 
apply: 


Momentum: 

DS = o 
Ds 

(y = Const. ) 

(3.3) 

Energy 

M = 0 

Ds 

(y = Const. ) 

(3.4) 

where 




A = Flow cross-sectional area in a streamtube 
C = Curvature of the streamline 
H = Stagnation enthaply 
n = Distance along the orthogonal 
r = Radial coordinate 
S = Entropy 

s = Distance along the streamline 
T = Static temperature 
V = Velocity 
y = Stream function 
p = Density 


The solution method is an extension of the conventional streamline curvature 
method. It is briefly described as follows: First a crude grid of streamlines 

and orthogonal lines is assumed ( see Figure 3.1). Second, the curvature of the 
streamlines at each of the grid points is evaluated. Third, the momentum equa- 
tion is integrated along a line normal to the streamlines to obtain velocity, 
and the continuity equation is integrated to determine the "correct" streamline 
positions (for the assumed streamline curvature field). These are indicated by 
the "x" in Figure 3.1. Fourth, an adjustment 6n, is computed by considering: 

(1) the difference between the computed and assumed streamline positions, and 

(2) the effect of the implied curvature modification in the integrated momentum 
equation. Finally, the streamlines are repositioned by the 6n values. 


Because the movement of any one grid point alters, through a change in curvature, 
the velocity at nearby points, it is desirable to account for these interrelating 
point adjustments simultaneously. The use of the simultaneous solution procedure 
employed in the STC program is not part of the classical streamline curvature 
method. In concept, the set of simultaneous equations for the normal stream- 
line adjustments are formulated from the finite difference equivalent of the 
following equations: 


a 2 ( 6n ) (1-M 2 ) a 2 ( 6 n ) = 

3Y 2 ( P V) 2 as 2 


(3.5) 
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where 


sn = Required streamline adjustment in the normal direction 
M = Mach number 

F = Driving (or error) function derived from the solution to the 
integral continuity and normal momentum equations. 


This equation is derived in Ref. 18 for the special case of isentropic two- 
dimensional flow. 



3 . 


4 . 


Integrate the cross- 
stream momentum equation 
and continuity equation 
to determine the "correct" 
streamline positions 


Solve the matrix equation 
for the an's and move grid 
points. 


Figure 3.1 STC Solution Technique 


THE TURBULENT BOUNDARY LAYER PROCEDURE 


The boundary layer procedure which has been coupled with the STC inviscid solu- 
tion is the method of Stratford and Beavers (SAB) which is described in detail 
in reference 19. In the SAB method, the integral boundary layer parameters, 
momentum thickness, displacement thickness, and boundary layer thickness are 
expressed in closed form as a function of Mach number, equivalent flat plate 
length, and Reynolds number based on equivalent flat plate length. 

Equivalent flat plate length is defined as the length over which a boundary layer 
growing on a flat plate at the given Mach number would acquire the same thickness 
as the real boundary layer at that given location. 
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For ax i symmetric or plane flow, this equivalent flat plate length is also de- 
fined in closed form as a function of Mach number and distance along the wall 
(with a radial coordinate in the axisymmetric case). 

Boundary layer separation is detected during the calculation by values of the 
Stratford separation parameter in excess of 0.5. 

On a single pass through the STC program, boundary layers are calculated (if 
required) at the end of the STC process, from substitution of derived flow para- 
meters into the closed form boundary layer relations described above. At this 
stage then, the STC inviscid flow field has not been merged with the computed 
boundary layers. 

The merging process is accomplished by restarting the program using the computed 
boundary layers from the first pass, and the flow near such boundaries is dis- 
placed accordingly. Several such restarts may be necessary before adequate con- 
vergence of the STC-SAB process is attained. 


OUTLINE OF PROGRAM STEPS 


The operations performed by the STC Program may be outlined as follows: 

(1) Define the flow regions and locate (approximately) the "primary" 
orthogonals and the streamlines which divide the internal and 
external flows. 

(2) Refine the grid as required by inserting additional streamlines and 
orthogonal lines between those already existing. 

(3) Compute the streamline angles and curvatures. 

(4) Compute the orthogonal line angles and move the grid points along the 
streamlines to obtain orthogonal ity. 

(5) Compute the velocities on the "far-field" boundary. 

(6) Adjust the flow rates in the exhaust streams, if any, to meet the 
calculated choking flow rate. 

(7) Integrate the momentum and continuity equations along each orthogonal 
[Eqs. (3.1) and (3.2)]. This step is also referred to as the "flow 
balance" calculation. 

(8) If the streamline positions are within a "rough tolerance", adjust 
the flow rate of variable flow channels to meet the "Kutta" condition 
at trailing edges. Step 7 is repeated with each iterative adjustment 
of the flow rates. 

(9) If the streamlines are within tolerance, return to step 2 for addi- 
tional grid refinement (unless grid refinement limits have already 
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been reached). Otherwise, continue to step 10. 

(10) Determine if the streamline positions are within final tolerance. If 
so, go to step 14. Otherwise continue to step 11. 

(11) Set up the matrix equation for the streamline correction, <5n. 

(12) Solve the matrix equation. 

(13) Modify the streamline position by sn, and return to step 3. 

(14) Calculate and print the output quantities; calculate boundary layers 
and adjust wake table at trailing edges; then return to step 1 for 
the next case, or for SAB restart. 

The first operation includes reading the input file for a description of the 
geometry and flow properties. The computer program has been written to have 
general capability for analyzing a great variety of configurations. The first 
step in the programmed logic is to develop a table of orthogonals or calcula- 
tion stations for the several flow regions. The regions are determined as 
illustrated in Figure 3.2 so the calculation can proceed from upstream to down- 
stream. The boundary of each region is defined as a primary orthogonal. As 
shown in Figure 3.2, the initial grid which is developed contains only the 
primary orthogonals and the double streamlines which separate the two illus- 
trated flow channels. 



Exterior 

Flow 

Channel 


Figure 3.2 Subdivision of Flow Field into Regions and Channels 

For a full description of successive steps in the computational procedure refer 
to references 17 and 18. 
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INTRODUCTION 


4 THE DUCT ACOUSTICS FINITE ELEMENT PROGRAM 


The finite element method described here was first implemented by Abrahamson 14 ** 5 . 
It is based on the usual perburbation of the basic fluid mechanics equations for 
small motions with mean flow parameters derived in the absence of fluctuating 
quantities and then substituted into the equations for the acoustic quantities. 
These equations are linearized by neglecting higher order terms and the problem 
is reduced to two dimensions by assuming mean swirl to be zero and assuming a 
harmonic solution in the circumferential direction. 

Originally, bilinear Serendipity isoparametric finite elements were used but the 
program has since been enhanced by the addition of subparametric Hermite ele- 
ments. The original Serendipity elements are conforming and of type C° (that is 
preserving continuity of variables across element boundaries). The Hermite 
elements are nonconforming C 1 elements preserving continuity of variables and 
approximate continuity of their first derivatives across boundaries. (A dis- 
cussion of the continuity characteristics of these elements is contained in the 
section on Element Derivation). 

Element matrices are formulated using a Galerkin procedure and assembled in a 
rectangular mesh for which the matrix map is specifically coded. Nodal vari- 
ables consist of the three complex components of acoustic velocity and (complex) 
acoustic pressure. In the case of Hermite elements, spatial derivatives of 
these physical parameters also become nodal variables. Nodal constants are the 
mean fluid parameters of density, velocity, and speed of sound. 

The resulting assembled global matrix is a block tridiagonal system in which the 
blocks are themselves block tridiagonal. The system is non-Hermitian and non- 
psoitive definite and is solved by forward and backward substitution following 
an L-U decomposition with pivoting restricted internally to major blocks. 

Since the solution of typical acoustic problems requires large solution grids 
and results in large demands on computer memory and central processor time, 
great attention has been paid to efficient memory management and efficient 
coding of computational operations. Use of the program has been simplified by 
the implementation of a fully automatic memory management scheme. 

For completeness, the theory described in the references is summarized here with 
a detailed description of program operation in its current form. 


DESCRIPTION OF THE MATHEMATICAL MODEL 


Derivation of Equations 


The nondimensional equations governing the motion of the fluid are: 
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I 


Conservation of Momentum 
+ 

p(|| + (V.v)vj = -vp + 1-v.i 


Conservation of Mass 


|| + v.(pV) = 0 


Conservation of Energy 

-p(iI + S.7f)- <r-l)(|| + V.5p) 

= — (— V.(icvf) + ( y- 1 ) 

Rg \ P i' / 


(4.1) 


(4.2) 


(4.3) 


Equation of State 

yp = ST (4.4) 


Where p, V, p, f , $, i are nondimensional density, velocity vector, pressure, 
temperature, dissipation function, and viscous stress tensor, respectively; 
y and k are gas specific heat ratio and thermal conductivity, and'R and P 
are Reynolds and Prandtl numbers. 

Equations (4.1) through (4.4) were made nondimensional for simplicity by using 
an arbitrary characteristic length; the ambient values of temperature, density, 
and speed of sound; and arbitrary characteristic values of viscosity and ther- 
mal conductivity. These equations govern both the mean flow and imposed acoustic 
motion within the duct. 

Although viscosity and dissipation were considered in deriving the mean flow, 
these influences were not considered in the acoustic equations. Thus, we write 
equations (4.1) and (4.2) in dimensional component form in cylindrical coordin- 
ates (ignoring viscosity) and assuming that all dependent variables are the sum 
of two motions, one the steady mean flow and the other a fluctuating acoustic 
motion, perturb these equations and subtract the unperturbed equations. 

Neglecting second order terms in the fluctuating quantities, assuming mean 
swirl to be zero, we get a set of first order equations in the fluctuating, or 
acoustic, quantities. 

The neglect of viscosity and heat conduction in the acoustic motion means that 
it is an adiabatic process. When such is true, the energy equation and 
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(4.5) 


equation of state may be manipulated to give 



Here c is the local speed of sound which is determined from the mean flow and p' 
and p' are acoustic pressure and density, respectively. 

Using equation (4.5) to eliminate p 1 and taking a harmonic solution in time (t) 
and angle (e) for the acoustic quantities, 

(u‘, v', w 1 , p ' ) = (u, v, w, p) e ut + m0 


where wand m are pure imaginary. The four equations governing the linearized 
acoustic motion in a nonuniform axisymmetric duct become 
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In these equations 0 , V, W and u, v, w are the mean and acoustic velocity com- 
ponents in the z, r, 9 directions, respectively. P and p represent mean pres- 
sure and density and p is acoustic pressure. 


Solution Strategy 


Two principal decisions in the adoption of a solution strategy were the use of 
a finite element algorithm and the choice of a uniform rectangular discretiza- 
tion mesh. 

A uniform rectangular discretization mesh does not imply that the grid points 
are to be equally spaced, or that the physical region should be rectangular. 
What it does imply is that there exists a simple mapping between nodes in a 
uniform rectangular mesh and nodes in an irregular nonuniformly discretized 
region. 
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The utility of this choice of a discretization mesh lies in the fact that it 
gives rise to an invariant global matrix map which may be explicitly embodied 
in a computer code for simplicity and efficiency. 

A uniform discretization mesh for a duct of constant cross-sectional area is 
shown in Figure 4.1(a). The matrix map generated for this discretization by 
a nodal numbering system such as that shown in the figure is developed below. 
Within each of the linear elements, interactions among nodal variables will 
only occur between those variables at the corners of the rectangle. Thus in 
element (2, 4), four sets of linear equations will link the four sets of nodal 
variables located at nodes 14, 15, 24, and 25. Consider the equation set rep- 
resenting variables at node 25. On assembly of the global matrix, this set 
will contain contributions from four adjacent elements namely (2, 4), (2, 5), 
(3, 4), (3, 5). That is, in the global matrix, the variables at node 25 will 
be linked explicitly with variables at nodes 14, 15, 16, 24, 26, 34, 35, 36. 
All other coefficients in the submatrix, representing the equation set for 
variables at node 25, will be zero. For the case of parameters per node, 
the matrix map generated by linear elements and the rectangular discretization 
numbered as shown in Figure 4.1(a), is given in Figure 4.2. 
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(b) Nonuniform Rectangular 
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(c) Nonuniform Rectangular Mapped by a Coordinate Transformation 


r 



Figure 4.1 Rectangular Discretization in a Plane. 
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Serial Node Number 
20 30 


1} Each submatrix Qhaa dimension "/x/ n 
where i = the number of parameters at a node. 

2) Matrix is Block Tridiagonal. 

3) Dimension of blocks = Number of Parameter! 
per row in serial numbering direction. 

4) Q)] denotes zero matrix. 



Figure 4.2 Matrix Map of the 2-D Rectangular 
Discretization of Figure 4.1 


The principal generalized characteristics evident from the map shown in Figure 
4.2 may be summarized as follows: 

■ Global matrix consists of three diagonal bands. 

■ Width of each band is "38-". 

■ Global matrix is block tridiagonal. 

■ Major blocks are themselves block tridiagonal. 

■ Order of major blocks is "n&", where "n" is equal to the number 
of nodes per row in the serial numbering direction. 

■ Order of minor blocks is "a". 

Note that an identical matrix map will be generated by the nonuniform 
discretizations in Figures 4.1(b) and 4.1(c). 
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Element Derivation 


The choice of a uniform rectangular discretization mesh imposes the requirement 
of a four-noded rectangular finite element as shown in Figure 4.3. 
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Figure 4.3 Rectangular Four-Noded Finite Element 


To allow this simple shape to conform to complicated boundaries, a linear trans- 
formation of the global coordinates is applied: 


r = (N lS N 2 , N 3 , N 4 ) 


z - (Nj , N 2 , N 3 , N 4 ) 



(4.10) 


where the (r., z. ), k = 1 to 4 are the global coordinates of the four corners. 
The functions N^are given by expressions of the form: 


N k = | (1 + ?5 k ) (1 + nn k ) 


(4.11) 


where values of U k , n k ) are fixed as (-1, -1), (1, -1), (1, 1), (-1, 1) 

1 to 4, respectively- K 

Within the deformed region, the mean fluid parameters p, U, V, c of equations 
(4.6) through (4.9) are also defined by relations of the type in equation (4.10). 
For example. 
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U = (Nj, N 2 , N 3 , ) / _2 



(N) (U> 


(4.12) 


Thus the element geometry and mean flow parameters are expressed in terms of 
nodal values and element local coordinates n and 5. 

Derivatives of the mean flow parameters are obtained as follows: 
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Thus for example. 


= (a ){U> 
3Z ' z' 


(4.14) 


The acoustic variables u, v, w, p within the finite element are given in terms of 
Hermitian functions of position and the values of the variables are as well as 
their derivatives at the nodes. For example, consider the variable u. 
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index "k": 
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H 1q (s) = 2s(s-l ) 2 

H u (s) = 2s 2 (s-1 ) 

where s = \ U+l) or s = i ( n+1 ) (4.18) 


Expressions of the type in equation (4.15) are inadequate for a finjte element 
representation, however, since the nodal variables (iiL) , ( 9U ) , ( 9 u ) are 

3? k % k Hrni'k 

referenced to local coordinates n and 5 . They may be expressed with respect to 
global coordinates as follows: 
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(4.20) 


(4.21) 


where the subscript "k" is used in every case to imply the value of the parameter 
at node k. 


Substituting equations (4.19) through (4.21) into equation (4.15) gives the final 
expression for the variables within the finite element: 
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(Note that in the above expressions, the M, ‘s only are functions of s and 
n.) k 
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Writing (4.22) as. 



(4.23) 


where the notation lu,} is used to represent u and its derivatives at ail four 
element nodes. We may define a matrix similar to (4.13): 




M - 1 1 (M) | 


(4.24) 


So that for example. 


9U 

3 Z 


(e z ) (u,} 


(4.25) 


Some discussion of the above expressions is helpful. In the rectangular element 
deformed with the n-C axis parallel to the r-z axis, for all values of 

k, and the contributions of nodal parameters (u,, r ) ana (u,, z ) are zero. 

Also, in this special case, continuity of derivatives u, r and u, 2 across element 
boundaries is preserved. This is not the case with any other element deformation 
and in general continuity of derivatives normal to element boundaries is not pre- 
served, although it always exists, of course, along element boundaries. 

Since solution of equations (4.6) through (4.9) requires continuity only of the 
variables themselves, it is thus possible to vary the number of nodal deriva- 
tives in (4.23) from 1 to 6. This option has been implemented in the computer 
program with the most commonly used options being 1, 3, or 4, yielding, respec- 
tively, a psuedol inear element, a nonconforming Hermite element, and a con- 
forming Hermite element for uniform area ducts. 

Applying the Galerkin version of the method of weighted residuals, we may inte- 
grate equations (4.6) through (4.9) over an element. This is performed numeri- 
cally using a Gaussian quadrature scheme in the simple undistorted local coor- 
dinate system. An element of area expressed in the local coordinate system is 
given by 


dA = | J | dndc 


(4.26) 


where |J| is the determinant of the Jacobian matrix. Equations (4.6) through 
(4.9) thus become: 
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Global Matrix Assembly and Boundary Conditions 


Assembly of a global matrix equation, that is, a matrix equation representing 
the entire system, from a set of finite element matrix equations is a basic pro- 
cedure in the finite element technique. Appropriate shifting of rows and col- 
umns are all that is required to add the local element matrix directly into the 
global matrix which is held in packed form (shown to the right in Figure 4.2) 
partially in direct access memory and partially on disk storage. 

Along the outer boundaries of the discretized region, several different kinds of 
boundary conditions may exist. The current procedure for incorporation of bound- 
ary conditions is the direct insertion into the global matrix of an equation 
representing the condition imposed on the boundary nodal variables. This pro- 
cedure is standard practice in finite element methodology, but has been shown 20- 
under some circumstances, to yield inferior results to the "natural" procedure 
for imposition of boundary conditions. The disadvantage with the natural pro- 
cedure is that boundary conditions are satisfied only in an integral sense over 
the boundary. Alternate procedures to both methods have been suggested, for ex- 
ample see the "characteristic correction" of Gottlieb et al 23 . These various 
techniques of imposing impedance boundary conditions are currently under review 
and alternate methods may be made available in the model. 

Also, under review is the procedure of applying noise "source" conditions by 
setting nodal pressure or velocity variables to a specified value. Here the 
specified value includes reflected as well as incident energy and separation of 
these requires additional analysis. The modal source and radiation conditions 
also presented in references 20 through 22 appear to offer a viable alternative 
to the current procedure and should soon be implemented in the model. 
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PROGRAM STRUCTURE 


Overview 


It is becoming generally accepted within the discipline of computer science that 
the best manner of presenting detailed documentation of a computer code is with- 
in the code itself. Accordingly, the acoustic finite element code is heavily 
commented with definitions of variable names, delineation of functional blocks 
and explanations of procedures being implemented. 

In contrast, the presentation given here describes the philosophy and overall 
structure of the program and is intended as an introduction for the general 
reader and programmer alike. 

The program is divided into four phases, namely, 

■ Data input, determination of memory structure, and memory management. 

■ Recovery and interpolation of aerodynamic data generated by the GE- 
STC Program. 

■ Preparation of an output file for receipt of results from the acous- 
tic finite element analysis. 

■ Execution of the acoustic finite element code. This involves two 
basic steps: 

1. Assembly and storage (on a disk file) of the global matrix 
equation. 

2. Reading global matrix (in portions from the disk file) followed 
by L-U decomposition and solution by forward and backward substi- 
tution. 

Memory Management and Input Data Structure 


Operation of the program is controlled by an executive which manages the sub- 
division of direct access memory into arrays for data storage. This procedure 
is very important in ensuring efficient use of memory and in permitting discret- 
ization grids of adequate size, since different phases of the program require 
different memory subdivisions. 

The procedure which is used is based on a principle embodied in the FORTRAN 
language whereby arrays are passed from a calling segment to a subroutine by 
a starting address only. 

Thus, using data input by the user, the executive calculates a memory subdivi- 
sion for each phase of the program. A Control Data Utility - Common Memory 
Manager - is directed to provide the required additional space, and array 
starting addresses are passed to subroutines for execution of each phase of the 
program. (This procedure is facilitated but not dependent upon Common Memory 
Manager. If this utility is not available on a particular computer installation, 
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then the required additional space must be provided explicitly by the user in 
"Blank Common".) 


Program input is described in detail in section 6, and consists of a data spec- 
ification subroutine and six user-supplied FORTRAN functions. These functions 
describe spatial variations in the density of the discretization, as well as 
admittance, source and radiation boundary conditions. 

Recovery and Interpolation of Aerodynamic Data 


Phase 2 of the program (embodied in subroutine AEROPRE) reads duct geometry and 
aerodynamic data from the random access file prepared by the GE-STC program. 

The aerodynamic data are, of course, defined on the mesh generated by the STC 
mean flow solution algorithm. Since this mesh contains partial streamlines and 
partial orthogonal lines, the logic for interpolation to the acoustic mesh needs 
careful implementation. Because the acoustic mesh has parallel radial lines, 
the method used consists simply of specifying the axial coordinate of one of 
these lines and interpolating along all streamlines defined there. At each 
acoustic axial station this yields a set of values unequally spaced in radial 
coordinate (r). Interpolation at specified points along these radial lines 
yields the acoustic mesh. 

Flow data output by the GE-STC program has not yet been merged with boundary 
layer flow profiles. This is performed subsequent to interpolation in AEROPRE 
by subroutine BNDLAYR and function BLPFILE using subroutine BLN0RM, from bound- 
ary thickness (6) data read from the GE-STC file. A one-seventh power boundary 
layer profile is applied to local Mach number so that the Mach number reaches 
.99 of its free stream value at a normal distance of 6 from the wall. 

Since mean flow temperature, pressure, density and the local speed of sound are 
all derivable from Mach number, ambient static pressure, and temperature, only 
Mach number and flow angle are interpolated and stored in array AERODAT. This 
array is retained in direct access memory throughout execution of the program. 
Subroutine AERO computes other mean flow parameters as required. 

Coordinates of the acoustic mesh are not stored but are computed as required 
from boundary coordinate data arrays by subroutine COORDS. 

Preparation of Output File 


Due to severe restrictions on direct access memory during the solution phase of 
the AFE program, the task of preliminary preparation of the output file for 
interactive graphics is performed somewhat out of sequence. The interactive 
graphics program requires the coordinates of the acoustic mesh and also Jacobian 
matrices for each finite element in order to perform interpolations. These two 
sets of parameters are computed, blocked together, and written out to a mass 
storage file by subroutine OUTPRE. 

Structure of the Finite Element Code 


The first step in the AFE code involves numerical integration of the finite ele- 
ment equations (4.27) through (4.30) to produce local element matrices which are 
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then assembled into the global matrix. Because of memory limitations, only 
the portion of the global matrix currently being assembled, is held in direct 
access memory. As soon as a global matrix block has been assembled with bound- 
ary conditions applied, it is written to a mass storage file. The procedure for 
performing the operations is described in more detail in figure 4.4. 

Defining the block structure of the global matric and its L-U decomposition fol- 
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The order of global matrix blocks on mass storage is then a f 1 b a 2 c f 

b n _ 1 a n _ 1 c n _ 2 f j b n a n c n _ 1 f p where the f ^ are vectors Trom the right-hand 

side of the matrix equation and vector f. corresponds to the i'th block row of 
the equation. 

The global matrix generated by the Galerkin finite element formulation is com- 
plex, non-Hermitian and nonpositive definite. Its L-U decomposition given above 
is readily obtainable by direct substitution: 


a k 


b k a k-i 

a k - b k c k-i 


for k = 2,n 


The system LUX = F now becomes LY = F and UX = Y, which may be solved by forward 
and backward substitution. 


*1 = f i 

■^k ~ ^k ^k ^k-i for k = 2,n 
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x k ■ a k l <*k - c k x k+l> for k ■ "- 1 ’ 1 

During assembly of the global matrix, the blocks are sparse and are packed in 
the manner described previously. During the solution phase "fill" is generated 
and an efficient minimum direct access memory configuration consists of three 
packed blocks, three block vectors and one expanded block. A description of 
the above operations using this memory configuration is given in Figure 4.5. 



Figure 4.4 Functional Diagram for Numerical Integration 
of Finite Element Equations 


28 






















5 INTERACTIVE GRAPHICS PROGRAMS 


INTRODUCTION 

Both the GE-STC mean flow program and the acoustic finite element program pro- 
duce large quantities of data, and in the case of both programs, it is necessary 
for the user to survey these data to satisfy himself as to the general "health" 
of the solutions over the entire domain. 

For this reason alone, good graphics capability is essential. A predetermined 
off-line graphics process is cumbersome and does not give the necessary flexi- 
bility for close examination of interesting features or anomolies that might be 
noticed during a cursory review. 

Interactive graphics, on the other hand, gives the user the capability to examine 
immediately an aspect of the solution behavior that might appear worthy of his 
attention. It is further possible to perform interactively, such functions as 
setting plot scales, adjusting the order of the interpolation between nodal 
values, or using the data to compute additional special parameters. 

The facility also exists within the two graphics programs for interactively 
specifying off-line creation of permanent high quality plots for selected 
cases. 

Flow Graphics Program 

The principal utility of the flow graphics program is to enable an ADAM user to 
ascertain rapidly whether the mean flow data file created by the GE-STC program 
is satisfactory. The program is relatively simple in construct, calculates no 
additional functions from input aerodynamic data and produces only three basic 
plots. Contouring, which in itself is a sophisticated process, is performed by 
a standard library package. 

The user may select contour intervals and plot scales and may examine selected 
areas of the flow field in detail. 

Plot types produced by the program represent contour maps in the r-z plane of 
three selected variables: 

■ Those streamlines which are continuous over the required intervals. 

■ Mach number contours at user-specified intervals. 

■ Velocity contours at user-specified intervals. 

Acoustic Graphics Program 

The acoustic graphics program gives a user the capability to plot amplitude or 
phase of any acoustic variable (three components of acoustic velocity or pres- 
sure) along any nodal line in the finite element mesh. Nodal acoustic data 
generated by the AFE program are stored on a file whose format is described in 


30 



Appendix 2. This file also contains data titles, important input parameters 
from the GE-STC and AFE programs, AFE-mesh coordinates and finite element 
Jacobian matrices. Selected portions of this information may be displayed on 
the interactive terminal screen to avoid user errors in data file identifica- 
tion. Between nodes, interpolation is performed using the following relation: 

4 

u(?,n) - + ^ 2 k 2 ^k^ ,n ^ + ^4k 3^k^ ,n ^ hir^i, 

k = l k 

+ ( J ik + J 3k 3 M k<^>) # k } 

taken from (4.22) with ^ set to zero for element boundaries. Interpolations 
are performed in local U,n) coordinates and transformed to global coordinates 
using the transformation, 

r = (N) Tr k > 

z = (N) {z k > 

where (N) = * [(1-5) (1-n), (1-0 (1+n), (1+0 (l+n),(l+c) (1-n)] 

The sequence of operations carried out in the program are listed below with pro- 
gram segments performing these operations shown in parentheses: 

1. Open random access file (MAIN) and read header records (RAHEADR) 

2. Set up dynamic storage areas (MAIN) 

3. Read nodal coordinates from random access file (RAINPT1 ) 

4. (User) specify whether to switch to radial or axial plotting mode 
(MAIN) 

5. For radial plots: 

■ Read in function values at nodes for selected axial plane 
(RCIN) 

■ Interpolate at user-specified radial locations for all avail- 
able acoustic variables and compute admittance if this plane is 
a boundary (RINTERP). 

■ Plot user-selected functions (RADPLOT) using NCAR graphics 
package (called from PLOTTR). 

6. For axial plots: 

■ Read in function values at nodes for selected radial nodal line 
(ACIN). 
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Interpolate at user-specified axial locations for all available 
acoustic variables, acoustic flux, and complex admittance if 
this nodal line is a boundary (ZINTERP). 


■ Plot user-selected functions (AXPLOT) using NCAR graphics pack- 
age (called from PLOTTR). 

The National Center for Atmospheric Research (NCAR) graphics package referred 
to above creates a "metacode file" which contains all completed plots. Special 
routines in the NCAR package combined with a NASA-Langley utility TRN and Tek- 
tronix graphics routines translate this metacode file to a display on the Tek- 
tronix terminal screen. 

Upon completion of an interactive plotting session, the user may direct TRN to 
access the NCAR "metacode file" and produce a "graphics command file" which is 
then accessed by an off-line batch stream job to produce permanent (CalComp) 
plots. 

Note: To ensure that direct access memory requirements are held below the 

limits allowable for interactive processing at NASA-Langley, this program 
performs dynamic memory management using the CDC utility "Common Memory Man- 
ager". Dependence on this utility may affect its portability. 
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6 USERS' GUIDE 


COMMENTS ON ADAM/NOS INTERFACE 


A more detailed version of the simplified block diagram of the ADAM system 
contained in Figure 2.1 is given in Figure 6.1. This figure shows in addition 
to the principal functional steps and data files, the CDC-NOS (Network Oper- 
ating System) interface. This interface consists of a set of procedures which 
greatly simplify the users task. Instructions on the use of these procedures 
are contained in this section for individual processes under the heading "Oper- 
ation Under CDC-NOS". 

Use of the system is strongly oriented towards interactive processes and in this 
context user ability to create, edit, and maintain NOS files is central to most 
operations. The principal tool for performance of these functions is the editor 
XEDIT. Manuals on NOS and XEDIT are available for users, from the NASA-Langley 
Analysis and Computation Division. 

Although functions, such as the creation of input data files, and the examina- 
tion of program output data are truly interactive, other processes are too con- 
suming of computer resources and/or take too long to be performed interactively. 
These include executing the GE-STC program, executing the AFE program, and the 
creation of high quality (CalComp) hardcopy plots. These processes are executed 
as "batch" jobs after being created by a user at an interactive terminal and are 
sent by him to the computer of his choice for off-line processing. 

Before executing ADAM NOS Procedures for the first time, a new user must get 
copies of these procedures (STCDECK, PROCPLT, AFEDECK, AFETDEK) saved under his 
user number and modified to include his own user and delivery information. 

RUNNING THE GE-STC PROGRAM 


Operation Under CDC-NOS 


Prior to execution of the GE-STC program, a "namelist" input data file should 
be prepared as described in the following sections using data sheets contained 
in Appendix 1. 

To execute the program the user then issues the following commands at an inter- 
active terminal 

GET, STCRUN 

STCRUN, DATIN = FI, DATOUT = F2, COMP = M. 

where FI = Name of user input data file 

F2 = Name which the user wishes to give to the output data file which 
will be created by the program. 

M = the letter name of the machine on which the user wishes the program 
to be executed. 
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Preliminary Comments on Input Data 


Primary items of importance in trouble-free use of the GE-STC program are: 

■ Specification of Smooth Input Geometry - Since the STC method uses 
boundary surface curvature to compute velocity and pressure gradi- 
ents, it is imperative that surface curvatures are piecewise con- 
tinuous. Smooth boundaries are most easily accomplished by using 
analytic functions with continuous second derivatives. 

■ Grid Refinement Criteria - The purpose of grid refinement criteria 
is to maximize refinement in areas of large flow gradients while 
minimizing additional refinement in areas of smooth flow. Exper- 
ience in program use will best assist the user in making optimum 
choices of refinement parameters when combined with the general 
guidelines in the following sections. 

■ Interation Tolerances - These control the logic for grid refine- 
ment and "flow balance" in the inner loop solution. Satisfactory 
operation usually occurs using the preset values. Instability in 
streamline adjustment may occur when small tolerances are combined 
with high grid refinement resolution in a region of large stream- 
line curvature. The remedy is to change tolerances, grid refine- 
ment or both. A good procedure is to keep iteration tolerances 
relatively large during grid development and then restart the last 
grid refinement with smaller tolerances. 

The STC analysis has been set up as a user-oriented program and most control 
variables have been preset to generally useful values. If numerical diffi- 
culties do occur, they will invariably be related to geometry errors or bad 
judgement in selecting grid refinement or iteration tolerances. 

INPUT DATA FOR THE GE-STC PROGRAM 


The following explanatory notes are intended to assist the user in completing 
input data sheets for the GE-STC program which are contained in Appendix 1. 
Although this program will cope with any consistent set of units, the Acoustic 
Finite Element Program expects values in the MKS system. Therefore, all input 
data should be specified in this system. (Note that values in parenthesis on 
the data sheets are default values to be used if no input is specified for 
particular parameters.) 

Input Sheet 0 

This consists of three lines of identification information in columns 2 through 
61 listing user name, user location, and problem identification, respectively. 

Input Sheet 1 


STC: Restart Control Parameter - Partially refined grids may be read from 

or saved on files by specifying a "T" in either column 14 or 24. In 
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the former case, it is necessary only to specify those input quan- 
tities which differ from those originally input. Note that the 
final run of a problem with boundary layers introduced into a fully 
converged STC solution is always a restart case. 

MACHO: Free stream Mach number 

2 

TSO, PSO: Ambient temperature (K), pressure (N/M ) 

RG, GAM: Gas constant J/(Kg.K) , ratio of specific heats 

RHL: Highlight radius (M) - A reference (or highlight) area (= tt(RHL) 2 , 

ax i symmetric or = RHL, planar) to be used for defining the mass 
flow for each channel (see AO sheet 3). 

RM: Maximum body radius (M) - Computed pressure drag forces are normal- 
ized by max. body area (= w(RM) 2 , axisymmetric, or = RM, planar). 

TTE: Body closure tolerance (M) - Maximum trailing edge thickness. 

AXI: Axisymmetric or planar (T or F). 

Grid Refinement Criteria 

The grid will be refined locally as required until the spacing of orthogonals 
and streamlines is less than the values given in the SGR and SGZ tables and the 
Mach number difference between any two points on a streamline or on orthogonal 
line is less than VMG1 and VMG2, respectively. As a guide, it is suggested that 
SGR be set to about twice the expected grid spacing. 

GR: Radial values where SGR is to be specified (M). 

SGR: Spacing of grid lines in radial direction (M). 

NGR: Number of values of GR and SGR. 

GZ: Axial values where SGZ is to be specified (M). 

SGZ: Spacing of grid lines in axial direction (M). 

NGZ: Number of values of GZ and SGZ. 

VMG1 , 

VMG2: Maximum Mach number increment between grid points in streamline 

or orthogonal directions, respectively. 

Iteration Tolerances 

The solution method consists of constructing a grid of streamlines and ortho- 
gonal lines. Starting with two streamlines per channel (one for each boundary) 
and an orthogonal passing through the first and last point of each boundary, 
the grid is automatically refined by successively dividing grid intervals in 
half. The solution accuracy and computational cost are directly related to the 
extent of grid refinement. Provided MAXIT, the maximum number of refinements, 
is not exceeded the grid will be refined locally as previously specified by the 
grid refinement criteria. 
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MAX IT : Maximum number of grid refinements. 

RHOC: Curvature damping - In some instances, it may be necessary to increase 

this parameter from its preset value of unity to a value of two or 
greater if an oscillation of MAXES2 is observed in the iteration 
history printout. 

NODENS: Number of streamline position refinements at constant density - At 

the default value of zero, the incompressible equations are used for 
the zeroth refinement only. At greater values (e.g. NODENS = 2) the 
inadvertant possibility of supersonic velocities occuring is pre- 
vented and is recommended for reliable starting. Where adjacent 
streams have large differences in total pressure (as in jet plume 
analysis), a value of -1 is recommended. 

Boundary Layer Parameter 

The following parameters define reference values for boundary layer computation. 
TREF: Reference temperature for viscosity solution (K) 

MUREF: Reference viscosity at TREF (J/Kg. K) 

SCON: Sutherland constand ( K ) 

Optional Input 

A large number of input quantities, in addition to those described above, may be 
specified. These items are not normally required for execution of the STC pro- 
gram and consist in general of controls for special program options and input 
data to modify preset or initialized constants or parameters. 

A selected few of these optional input parameters are described in the addendum 
to Appendix 1. For a full list, the user is directed to Reference 17. 

Input Sheet 2 


Each input sheet 2 specifies data along one particular boundary. When preparing 
an aerodynamic data file for use by the AFE program, only two geometric bound- 
aries are permitted, namely inner and outer duct walls. To get accurate esti- 
mates of flow for real engine nacelles, however, it is often necessary to per- 
form preliminary STC program runs with several flow channels and multiple boun- 
daries. For example. Figure 6.2 shows a case with three flow channels. Upon 
completion of the preliminary analysis, the positions of stagnation streamlines 
will be available. These may then be used as artificial boundaries for a single 
channel analysis for use by the AFE program. 


37 


Far Field Boundary 

n 



Figure 6.2: Illustration of a Multiple Channel Configuration 


Each channel must be bounded, at least in part, by an upper boundary and a lower 
boundary. It is also possible for a channel geometric boundary contour to be 
composed of several physical boundaries. In this case, an input sheet must be 
completed for each physical boundary and end point coordinates must be identi- 
cal. The AFE program will also recognize this situation and will perform cor- 
rectly, linking physical boundaries to form a single geometric contour. Certain 
channel and boundary names are reserved by the program for particular uses as 
listed below. 

Reserved Channel Names: 

EXT - External flow channel 
INT - Inlet flow capture channel 

Reserved Boundary Names: 

FF Far field boundary. This name causes the numerical flow 

solution to be matched to a small pertubation analytical 
solution in the region from the far field boundary to in- 
finity. 

PRES II Up to two pressure boundaries are allowed permitting static 
PRES 2 J pressure along a boundary to be specified by the user. 

FREE1 I A free boundary is a constant pressure boundary which is 

FREE2 ( downstream of a fixed boundary and has the same pressure as 

the last point on the fixed boundary. 

Items listed on input sheet 2 are generally self explanatory. Boundary co- 
ordinates at sharp corners must be listed twice, one time for each angle which 
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exists at that point. In each interval, a locally rotated cubic polynomical 
is fitted by the program. Generally, it is extremely important that contours 
are smooth and points consistent with specified angles. 

All points must be listed in the streamwise direction. For an inlet lip 
starting at the highlight, internal and external surfaces are listed separ- 
ately with ANGD at the highlight equal to +90° for the external surface and 
-90° for the internal surface. 

The format for specifying pressure boundaries or free boundaries is contained 
in input sheet 2a. 

Boundary layer calculations on a particular surface, are normally specified 
using the input given on sheet 2. In some cases, it is desireable to intro- 
duce boundary layers into a fully converged STC inviscid solution. This may 
be conveniently done on restart using the procedure described in input sheet 
3a. 

Input Sheet 3 


Data on input sheet 3 supplies entrance flow properties for each channel. 

Total pressure and total temperature may be specified using one of the follow- 
ing procedures: 

■ Specify TTO and PTO if the stagnation properties are known. 

■ Specify MACHO, TSO and PSO if static properties and Mach number 
are known. If only MACHO is supplied, values from sheet 1 will 
be used for TSO and PSO. 

■ If no values are input, then free stream values from sheet 1 are 
used. 

AO is the channel area used for determining channel flow rate. It should be 
input as a fraction normalized by the highlight area calculated from RHL on 
sheet 1. 

If for any channel input data are not supplied, the reference properties on 
sheet 1 will be used with frontal area computed at the entrance station. This 
option is suggested for an external stream. 

Channel flow rates will normally be adjusted to meet trailing edge or choked 
flow conditions. If flow rate is not to be varied for this channel, set VARY 
F. 

An example of a typical input data file for the GE-STC program is contained in 
Appendix 1 following input data sheet blanks. 

I nterpretation of Program Output 

Output from the STC program may be divided logically into the following eight 



sections 


1. Input data file. 

2. Input and calculated boundary coordinates and angles. 

3. Printer plot of initial XII, XI2 grid map. 

4. Solution history. 

5. General input and output data. 

6. Flow field data along orthogonal lines. 

7. Calculated flow data along field boundaries and final channel 
momentum balances. 

8. Boundary layer data. 

Only items 6, 7, and 8 will be described in detail, as the other items are 
best left to the reference. However, a few notes are in order. Problems 
arising in the use of this program can usually be attributed to incorrect 
specification of the initial or boundary conditions. A careful inspection of 
the calculated boundary quantities of section 2 will usually show a point or 
region of high curvature, which should be examined for correctness. Also, 
in the solution history, the behavior of MAXES2 relative to LIMES2 indicates 
convergence by a decreasing absolute value. Any erratic behavior in MAXES2 
indicates convergence problems and the user is referred to notes in the adden- 
dum to Appendix 1 and to reference 17. 

6 Flow Field Data Along Orthogonal Lines 

A tabulation of column headings in the section of data output follows: 

XII - Orthogonal ( station) numbering coordinate 

XI2 - Streamline numbering coordinate 

STRM FNCT - Dimensionless stream function (fraction of flow in channel). 


X, Z - Streamwise coordinate 

Y. R - Transverse or cross-stream coordinate 


PHI - 

Flow angle (degrees) 

CURV - 

Streamline curvature 

PS/PO - 

Ratio of local static 

PS/PT - 

Ratio of local static 

TS/TT - 

Ratio of local static 

CP - 

Pressure coefficient 

MACH - 

Mach number 

Area - 

Flow area 


pressure to reference static pressure 
pressure to total pressure 
temperature to total temperature 
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PT/PTO - Channel total pressure ratio 

1 _ Calculated Flow Data along Field Boundaries and Final Channel Momentum 
~ Balances 


S1W - 

Distance along boundary 

xw,zw - 

Axial coordinate 

YW,ZW - 

Vertical or radial coordinate 

ANGW - 

Flow angle or surface angle 

CURVW - 

Streamline or surface curvature 


Momentum Balance - The difference between the entering and leaving momentum 

may be interpreted as an indication of the accuracy of 
the analysis. 


8 Boundary Layer Data 

XW - 

Axial coordinate 

THETA - 

Momentum thickness (e) 

DSTAR - 

Displacement thickness (s*) 

DELTA - 

Boundary layer thickness («) 

REX - 

Local Reynolds number 

CAPX - 

Equivalent flat plate distance along surface 

CF - 

Skin friction coefficient (C f ) 

SW - 

Distance along surface (s) 

DSTR - 

Smoothed displacement thickness (6*) 

DDSTR - 

ds*/ds 

SEP - 

Separation flag, = SEP if separation occurs 

F - 

Stratford separation parameter 

USING THE 

FLOW GRAPHICS PROGRAM 

Operation 

under NOS 


After a random access mass storage file containing aerodynamic data from the 
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GE-STC program has been generated, it may be examined interactively by issuing 
the following commands from a Tektronix terminal: 

GET, TAPE 11 = (filename of GE-STC output file) 

GET, PROCPLT. 

BEGIN, STRPLOT, PROCPLT. 

The procedure FEPLOT executes the GE-STC data graphics program and subsequently 
under user direction, may proceed to execute a post processor program to allow 
the user to obtain CalComp permanent plots via a batch job submittal. 

The postprocessor program which creates the "Plot Vector File" at NASA-Langley 
is called TRN. The use of TRN for permanent plots from the flow graphics pro- 
cedure is identical to its use by the acoustic graphics procedure and is de- 
scribed under that heading. 

Notes for an Interactive Session 

Throughout the interactive session, all user input is prompted by directions 
appearing on the screen. 

For all user input, a "Y" indicates "yes" and "N" indicates "no". A blank or 
any other character is assumed to be "no". 

Queries are self explanatory and when user input data are required they are 
checked to ensure that they fall within the valid range before proceeding. 

The three plots are produced in sequence so that in order to repeat a particu- 
lar plot, it is necessary to repeat the entire sequence. 

When completing an interactive session, the user may enter the postprocessor 
program TRN. No prompting instructions are issued by TRN and the user is obli- 
ged to know the required command format. 


RUNNING THE ACOUSTIC FINITE ELEMENT PROGRAM 


Operation Under CDC-NOS 

A convenient procedure has been devised for running the AFE program under the 
CDC-NOS operating system: 

Step 1 - Creation of a user-functions input file as described in the 
following section. 

Step 2 - Creation of a geometry and aerodynamic data file using the 
GE-STC program. 

Step 3 - At an interactive terminal, the user issues the following 
instructions: 


GET,AFERUN. 



AFERUN, FUNC = Fl, AERO = F2, OUT = F3, COMP = M. 


where Fl = User-functions file name. 

F2 = The name of the output file which was created by the GE-STC Program. 

F3 = The name which the user wishes to give to the output file which will 

be created by this run of the acoustic finite element program. 

M = The letter name of the machine where the user wishes the program to 

be executed. 

If the user desires a listing of the program or has made changes to the code, 
then the following instructions should be issued: 

GET,AFETRUN. 

AFETRUN,PROG = NAME, FUNC = F2, OUT = F3, COMP = M. 

where "NAME" is the name of the new program source file. 

The ADAM NOS procedures described above use the following files: 

AFEDECK - The job deck which will be sent to the computer of the users' 
choice to run his job. 

AFERUN - The interactive procedure to send the job to a batch processing 
stream. This procedure also performs a syntactic check on the 
user functions file, and checks to see whether the GE-STC file 
exists. 

AFETDEK - The job deck to test run a new version of the AFE Program. 

AFETRUN - The interactive procedure to send a test run to a batch processing 

stream. In addition to the actions performed by AFERUN, a syntac- 
tic check on the new version of the AFE program is also performed. 

AFE 12 - The program source file. 

PBIN - The compiles relocatable binary file from the source file, AFE12. 
Preparation of the User-Function File 

User input to the AFE program consists of a FORTRAN subroutine and five FORTRAN 
functions. 

Subroutine DATASUB - Specification of input variables. 

Function ZPDFN - Axial points distribution 

RPFN - Radial points distribution 

Outer wall acoustic admittance 


ADMITO - 
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oooooooooooooooooooooooooooooooo ooo 


ADMITI - 


Inner wall acoustic admittance 


ADMIT!" - Admittance of duct termination 

PSOURCE - Pressure distribution over source plane 

The FORTRAN subroutine DATASUB is simply a vehicle for transferring data into 
the program. This could have been equally specified as a data file but has 
been cast as a subroutine for consistency in format with other program input. 
In this way, a single file may be readily used for all input data. 

The format of the data subroutine is shown in the following example: 

SUBROUTINE DATASUB (M,NDUCT,ZLOCS,ZLOCR, FREQ, MSP IN, FLOW 
1 ,AMPHAS, OUTFILE,BCWEAK,NNDERIV, IWRITE) 

This subroutine contains input variables for program AFE12. 

LOGICAL FLOW,AMPHAS,OUTFILE,BCWEAK 

Set Program Input Variables: 

M = No. of elements in radial direction 

NDUCT = No. of axial elements in duct. 

ZLOCS = Z-coordinate of source in GE-STC reference frame 

ZLOCR = Z-coordinate of radiation boundary in GE-STC reference frame 
FREQ = Source frequency in hertz 

MSPIN = Circumferential (spinning) mode number 
= 0 for axisymmetric case. 

FLOW = .FALSE, for zero mean flow 

AMPHAS = Output format indicator 

= .TRUE, for output of nodal variables in amplitude and phase 
= .FALSE, for output of nodal variables in real and imaginary 

Outfile = .TRUE, for creation of output file containing analysis 
results for use by interactive graphics program. 

BCWEAK = .TRUE, for application of derivative wall boundary conditions 
by subroutine LINK06. 

= .FALSE, for successive applications of derivative wall boundary 
conditions by subroutine LINK02. 

NNDERIV = Number of nodal derivatives of each physical parameter 
= 1 for linear elements 
= 3 or 4 for Hermite elements 
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C IWRITE = 0 for no printed output 
C 1 for nodal variables only 

C 2 for nodal variables + interpolated values in radial dirn. 

C DATA TITLE : PROGRAM TEST FOR GEDATA2 

M = 11 
NDUCT = 50 
ZLOCS =3.0 
ZLOCR = 0.01 
FREQ = 512.0 
MSPIN = 0 
FLOW = .TRUE. 

AMPHAS = .TRUE. 

OUTFILE = .TRUE. 

BCWEAK = .FALSE. 

NNDERIV = 3 
IWRITE = 2 
RETURN 
END 

Where lines beginning with a "C" are nonexecutable comments, explaining the 
function of input variables. 

Variables "M" and "NDUCT" define the overall size of the mesh. The magnitude 
of M is limited by the size of direct access memory, whereas NDUCT is unlimited. 
Current limits to M on a Cyber 175 computer with just over 100k words of direct 
access memory is 46 linear elements or 12 Hermite elements. 

Spacing of grid points in the axial and radial directions is controlled by 
points distribution functions ZPDFN and RPFN, respectively. These functions 
define curves of normalized coordinate distance versus point number. 

For example, the following function specifies a uniformly spaced axial mesh by 
giving a linear relationship between ZPDFN ( =Normalized axial coordinate) and 
axial point number J, where J varies from 1 to N+l. (Note that for N finite 
elements there exist N + 1 points). 

Function ZPDFN (J,N) 

C 

C ZPDFN: Normalized axial (Z) points distribution function 

C J = Axial point number 

C N = Number of axial elements (= No. of points - 1) 

C ZPDFN = Normalized axial coordinate corresponding to point J 
C 

C This is a linear points profile 
C 

ZPDFN = ( J— 1.0) /N 

RETURN 

END 

More complicated distributions are frequently necessary. For example, analysis 
of a duct with flow containing a constriction in the central area of the analy- 


45 


sis region requires a finer mesh spacing in the constriction region since acous- 
tic wavelength decrease here for sound propagation against the mean flow. 

The following function (illustrated in Figure 6.3 for N = 50) was used for this 
problem. 

Function ZPDFN(J,N) 

C 

C ZPDFN: Normalized axial (Z) points distribution function 

C J = Axial point number 

C N = Number of axial elements (= No. of points - 1) 

C ZPDFN = Normalized axial coordinate corresponding to Point J 
C This is a non-uniform profile with points clustered about the center 
C of the Z-coordinate range. 

DATA PI/3.14159265/ 

JM1 = J-l 

IF (JM1.GT.N/2.) GO TO 10 
ZPDFN = 0.5*SIN (JM1*PI/N) 

RETURN 

ZPDFN = 0.5*(2.+SIN(- (N-JM1 )*PI/N) ) 

RETURN 

END 



6.3 Nonuniform Axial Points Distribution Function for 
Duct with Central Constriction 


In the radial direction the presence of a boundary layer on the outer wall with 
large flow gradients requires a closer spacing of mesh points. A points distri- 
bution function specified explicitly for such a case with M=ll follows: 
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Function RPFN(I,M) 


C 

C RPFN - Normalized radial points distribution function 
C I = Radial point number (input) 

C M = Number of radial elements (input) (= No. Points - 1) 

C RPFN = Normalized radial coordinate corresponding to Point I 
C 

C This is an example of a radial points profile with data specifically 

C set for 11 radial elements. 

C 

DIMENSION DAT (12) 

DATA DAT/0.0,.01,.1,.3,.5,.7,.9,.99,.999,.9999,.99999,1./ 

IF (M.NE.ll) STOP 22 
RPFN = DAT (I) 

RETURN 

END 

In this profile, a single point is also placed close to the axis, although no 
boundary layer exists here, since this is often a region of large acoustic 
gradients. 

A more general profile accomplishing a similar end to the previous case follows: 

Function RPFN(1,M) 

C 

C RPFN = Normalized radial points distribution function 
C I = Radial point number (input) 

C M = Number of radial elements (input) (=No. Points - 1) 

C RPFN = Normalized radial coordinate corresponding to Point I 
DIMENSION DI(3), 00(8) 

DATA DI/0.0, .005,-015/ 

1, DO/. 96,. 98,. 993,. 998,. 9995,. 99995,. 999995, 1.0/ 

C 

C This DISTRIB FN. is composed of three segments with inner and outer 

C segments specified explicitly and the middle segment being a linear 

C division among remaining elements. 

C 

5 NI = 3 

NO = 8 

C NI, NO = Number of nodes in inner and outer distributions, respectively. 

IF (I.GT.NI) GO TO 10 
RPFN = DI(I) 

RETURN 

10 IF ( I.LT.M-NO+2) GO TO 20 

RPFN = DO ( I-M+NO-1 ) 

RETURN 

C NC = Number of elements in central distribution 
20 NC = M-NI-NO+2 

RPFN = DI(NI) + (D0(1)-DI(NI))*(I-NI)/NC 

RETURN 

END 
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Complex nondimensional specific acoustic admittances on boundaries of the analy- 
sis region are specified by the three functions ADMITO, ADMITI, AND ADMITT. 

Here ADMITO refers to the outer wall boundary, ADMITI to the inner wall bound- 
ary for a coaxial duct, and ADMITT to the termination boundary. The format of 
these three functions is identical; all being dependent on a single nondimen- 
sional spatial variable. Note that the analysis employs the e + J w t convention 
so that when comparing admittances with analyses using the e-jwt convention, 
the sign of the imaginary part should be reversed. Note also that standard 
values of p and c are used in converting to dimensional admittances. These 
values are specified at the beginning of the program in MKS units. 

A "hard wall" example of an admittance function follows: 

COMPLEX FUNCTION ADMITI (Z ) 

C 

C ADMITI = Specific acoustic admittance of inner walls 
C Z = Nondimensional axial distance (Range is from zero to one). 

ADMITI = CMPLX(0. 0,0.0) 

RETURN 

END 

An example of a "plane wave" radiation admittance (= ^) follows: 

COMPLEX FUNCTION ADMITT(R) 

C 

C ADMITT = Specific acoustic admittance of duct termination 
C R = Nondimensional annular radial coordinate (= 0 at inner wall) 

C 

ADMITT = C0MPLX( 1 .0,0.0) 

RETURN 

END 

The function PS0URCE allows the user to specify a complex pressure distribution 
over a radial plane at one end of the analysis region. 


e.g. 

COMPLEX FUNCTION PS0URCE(R,MSPIN) 

C 

C Function calculates pressure distribution at duct source 

C R = Nondimensional annular radial coordinate (= 0 at inner wall) 

C (= 1 at outer wal 1 ) 

PS0URCE = CMPLXd. 0,0.0) 

RETURN 

END 

Note that this pressure distribution includes incident acoustic energy, as well 
as acoustic energy reflected back from the duct itself. 

USING THE ACOUSTIC GRAPHICS PROGRAM 
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Operation Under NOS 

After a random access mass storage file containing acoustic data from the AFE 
program has been generated, it may be examined interactively by issuing the 
following commands from a Tektronix terminal: 

GET,TAPE1 1 = (filename of AFE output file) 

GET, PROCPLT. 

BEGIN, FEPLOT, PROCPLT. 

The procedure FEPLOT executes the AFE data graphics program and subsequently 
under user direction, may proceed to execute a postprocessor program to allow 
the user to obtain CalComp permanent plots via a batch job submittal. 

The postprocessor program which creates the "Plot Vector File" at NASA Langley 
is called TRN. Its input is prompted by "TRN>" on the screen. A document 
"Introduction to TRN" may be obtained interactively as follows: 

GET, TRNDOC/UN = 973650N. 

Related files in this section of the ADAM system are: 

FEPLFOR - Interactive graphics program source file. 

FEPLTAB - Compiled binary executable file of FEPLFOR. 

Guide to a Typical Interactive Session 

Throughout the interactive session, all user input is prompted by directions 
appearing on the screen. Since these directions are self-explanatory, these 
notes are brief. 

The first action of the program is to write titles and a list of data-identify- 
ing parameters on the screen. These are obtained from the AFE data file and 
relate to both GE-STC program input parameters and AFE program input parameters. 

The user is then given the option of verifying AFE-mesh nodal coordinates at 
boundaries. 

Subsequently, the graphics portion of the program begins and the user is able 
to enter either the radial or axial plot mode. 

In either mode, nodal coordinates at the specified nodal index are displayed 
and parameter plots along the entire nodal line may be made using nodal values 
only, or nodal values with a specified number of interpolated values between 
nodes. The capability for plotting only a segment of the nodal line also 
exists and appropriate end points may be specified. 

The number of available plot options depends on position and orientation. The 
following options are always available in both amplitude and phase: 

1 - Axial velocity 

2 - Radial velocity 
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3 - Tangential velocity 

4 - Pressure 

In addition, in the axial direction plot mode a fifth option, acoustic flux 
integrated over that radial plane is available, and is the same regardless 
of the radial node specified. 

Also at any boundary, acoustic admittance may be computed and its real and 
imaginary parts plotted. 

A "no plot" option allows the user to progress to another nodal index station 
or to change directional mode. 

With every plot, scales are set automatically but the user may change them if 
the resultant plot is not to his liking or to demonstrate a particular facet 
of the data. 

Generally, the following instructions apply to all user input: 

■ "Y" indicates yes 

■ "N" indicates no 

■ "?" displayed on the screen with no prompting directives allows 

the user time to make a hardcopy of the screen contents. To 
proceed to the next step, enter "carriage return". 

At the conclusion of an interactive session, a directory of all plot frames 
created during the session are displayed. 

If permanent CalComp plots are required, the user may enter the postprocessor 
program TRN. 

Within TRN, no prompting directives are issued and the user is required to know 
the format of TRN instructions. 

For simple use, only three commands are required: 

P - Plot 

5 - Skip 
G - Go 

Q - Quit 

P and S are followed by a blank space and a number indicating the number of 
frames the user requires plotted or skipped. 

Input commands are separated by blanks, commas or the end of a line. 

Thus for example, to plot frames 2, 3, 4 and 7, 8 out of a ten-frame direc- 
tory, the following instructions would be issued in response to the TRN> 
prompt: 
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S 1, P 3, S 2, P 2, G 
Then response to the next TRN> prompt: 

Q 

Procedure FEPLOT then submits an off-line batch job for creation of the re- 
quired plots. 

Examples of plots created using this process are shown in Section 7. 
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7 EXAMPLES AND DISCUSSION 


MODEL CONFIGURATION 

The configuration selected for testing and demonstration of the ADAM system 
is relatively simple. It consists of a narrow variable area duct with hard 
walls, approximately three meters long with uniform cylindrical sections at 
each end tapering gently to a centrally located constriction (see Figure 7.1). 
The radius of cylindrical sections is 0.15 meter and the area ratio of these 
sections to the constriction throat is =2. 


Radiation 

Plane 

c 


Constriction 

Throat- 


Source 

Plane 


Flow 


c 


1 

0.3m 

T 


Meters 

Figure 7.1 Test Case Geometry: Duct With Central Constriction 

This configuration has already been analyzed for the "no flow" case by Silcox 
and Lester-^. They found that results from the AFE program agreed well with 
measured values over a range of modal acoustic source distributions and fre- 
quencies. 

The purpose of the following examples is to illustrate the format of the ADAM 
system graphical output and to show that the system produces anticipated re- 
sults when applied to a simple configuration. 

APPROXIMATE ENERGY ANALYSIS 

In order to provide an independent assessment of sound propagation character- 
istics in the duct configuration described above, the following approximate 
analysis is derived from conservation of energy principles. 

Consider the acoustic intensity in a variable area duct with quasi-irrotational 
flow. Morfeys expression for the intensity vector (see references 24 and 25) 
is: 

N = pV + p($ 0 .$)$ + J- 2 $ 0 p 2 + -1 $,(0)p (7.1) 

pC c z 

where N and V are the acoustic intensity and velocity vectors, respectively, and 
$ o is the mean velocity vector. 
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Expanding this expression in terms of its components, taking the time average 
and assuming mean swirl to be zero, we get the following expression for the 
time average of the axial acoustic intensity <N >: 


<N z > = |Re {(1+M*) pu* + M z M r pv* + -1^- M z pp* 
+ pcM z uu* + pcM r uv*> 


(7.2) 


where M r ,M z are the radial and axial Mach numbers, "Re" indicates "real part 
of" and indicates complex-conjugate. 

This is also the expression presently contained in the AFE program and is in- 
tegrated over a radial plane for computing acoustic flux. 

Suppose now that the duct contains a plane wave propagating axially against 
the mean flow and that minimal reflections are caused by the constriction 
and duct termination. Since the constriction taper is very gradual, the 
radial Mach number ( M ) is a second order quantity as is radial acoustic vel- 
ocity (v). 

Under these conditions, acoustic pressure (p) and acoustic axial velocity (u) 
are in phase and are simply related by the expression: 



Equation (7.2) then becomes: 
<N Z > = It ( 1+M z )p} 2 


Acoustic flux (F) is then simply the product of this expression and the cross- 
sectional area. Since the duct has hard walls and no other energy absorbing 
mechanisms have been introduced, net flux is conserved with distance down the 
duct. Setting the flux in the constant area section of the duct equal to the 
flux at the throat gives an expression for the ratio of the pressures at these 
two locations: 


e _t n + yr c (7 

P c ‘ (1 + M Zt >r t 

where subscript "c" refers to constant area, subscript "t" refers to throat, 
and r is the radius. Similarly, we find that the expression for acoustic 
axial velocity ratio (u./u ) is identical to that given for the pressure in 
equation (7.3). 
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STC ANALYSES 


Two STC analyses with different flow Mach numbers were performed. Input data 
to the STC program for the higher Mach number case are contained in Appendix 1. 

Plots of streamlines, Mach number contours and velocity contours for this case 
obtained from the flow graphics program are shown in Figures 7.2 through 7.4, 
respectively*. Similar plots for the high Mach number case are shown in Fig- 
ures 7.5 through 7.7. Dashed lines (with long dashes) in all flow graphics 
figures indicate boundary layer thickness while lines constructed of short 
dashes, on streamline plots only, represent boundary layer displacement thick- 
ness. 


CORRELATION OF APPROXIMATE ENERGY ANALYSIS WITH RESULTS OF 
AFE ANALYSES 

Using equation (7.3) with the results of the STC analyses, we may determine 
acoustic pressure ratios between cylindrical duct sections and constriction 
throat for the two flow cases described above. These pressure ratios, to- 
gether with that for a zero Mach number, are shown in Table 7.1. 



Approx. Energy 
Analysis 

AFE Program 
Results 

Case 

Number 

c 


_P t *4 

pT (or TT> 
^C c 

— Fj u t 

Pc u c 

1 

0 

0 

1.43 

1.45 1.41 

2 

-.141 

-.295 

1.74 

1.70 1.80 

3 

-.27 

-.69 

3.4 

2.9 3.6 


Table 7.1: Comparison of Acoustic Pressure 

and Velocity Ratios from Approx- 
imate Analysis with AFE Analyses 


The Mach numbers shown in Table 7.1 represent average Mach numbers across a 
radial plane. With regard to these values, it is interesting to note that 
although the area ratio of constant duct section to constriction throat is 
equal to 2, the ratios of the Mach numbers only approximate to this value 
(i.e. 2.1 and 2.6). This is the result of compressibility in the flow pro- 
gram and its effects become more pronounced at higher Mach numbers. 

*Figures 7.2 through 7.27 follow the text in Section 7. 
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Correlation for No Flow 


From the preceeding approximate analysis and from the results presented in the 
above table, we may thus expect that a plane wave propagating down the duct in 
the absence of mean flow would exhibit equal uniform pressure amplitudes in the 
constant area sections increasing to 1.43 times their magnitude at the throat of 
the constriction. Also, the axial component of acoustic velocity would be ex- 
pected to display identical characteristics. 

Results from an AFE analysis using linear elements with a mesh consisting of 20 
radial elements by 400 axial elements are shown in Figures 7.8 through 7.16. 
Source pressure distribution is uniform and frequency is 512Hz, well below the 
cut-on frequency for higher order radial modes. 

Figures 7.8 and 7.9 show the axial dependence of acoustic pressure, amplitude 
and phase, along a line located approximately midway between duct centerline 
and outer wall. 

As expected, pressure in Figure 7.8 is approximately constant in constant area 
sections of the duct and peaking near the constriction throat with an amplitude 
close to the value predicted by the approximate analysis (see Table 7.1). The 
parallel straight phase lines in Figure 7.9 are characteristic of plane wave 
propagation. Clearly, the approximate analysis is not quite correct in assuming 
no reflections from duct constriction and termination as shown by the standing 
wave pattern superimposed on the basic pressure characteristic in Figure 7.8. 

Axial dependence of the amplitudes of axial and radial velocity components 
along the same line is shown in Figures 7.10 and 7.11, respectively. Axial 
velocity displays a similar trend to pressure (Figure 7.8), while radial 
velocity is nonzero only in nonuniform sections of the duct. Interestingly, 
it tends to drop quickly towards zero as the duct walls briefly become par- 
allel at the constriction throat. 

Other assumptions of the approximate analysis are verified by the results in 
Figures 7.12 through 7.18. Acoustic flux integrated over radial cross-sections 
is shown to be conserved (Figure 7.12), a plane wave evidenced by constant 
radial pressure profiles is seen to be propagating in uniform sections of the 
duct (Figures 7.13 and 7.14), and acoustic admittance of the duct outer wall 
is zero (Figures 7.15 and 7.16). 

Correlation With Flow 


For the case of sound propagation along the duct in the opposite direction to 
the mean flow the approximate analysis shows that the relative heights of the 
pressure, and axial velocity peaks at the duct constriction throat increase 
with flow velocity (see Table 7.1). 

AFE analyses using Hermite elements with first order derivatives only, on a 
mesh of 11 radial elements by 50 axial elements for the two nonzero flow 
cases in Table 7.1 were performed using the axial points distribution func- 
tion in Figure 6.3 with the explicit radial points distribution function 
and the DATASUB example also shown in section 6. Results are presented in 
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Figures 7.17 through 7.26. 


Figures 7.17 and 7.21 show axial dependence of acoustic pressure amplitude along 
a line approximately midway between duct centerline and outer wall for the two 
nonzero axial flow velocities in Table 7.1, respectively. Ratios of throat to 
uniform section pressures are seen to increase with flow velocity to a peak 
agreeing closely with values predicted by the approximate model (see Table 7.1). 
Axial velocity components along the same line increase in a similar manner as 
shown in Figures 7.19 and 7.23. 

Comparing Figures 7.9, 7.18, and 7.22, axial dependence of acoustic pressure 
phase shows the expected effect of increased mean flow velocity by an increased 
rate of change with distance. This axial rate of change of phase with distance 
corresponds to the propagation constant (k) in a one-dimensional model. Based 
on a one-dimensional model, it would be expected that k would have the follow- 
ing dependence on Mach number (M): 



1-M 


where k Q is the axial propagation constant for zero flow. 

The effect of this nonlinear dependence on Mach number is seen by the magnified 
effect of the constriction in the higher Mach number case (Figure 7.22) com- 
pared to the lower Mach number case (Figure 7.18). 

Acoustic flux is again conserved in both cases as shown in Figures 7.20 and 
7.24. 


CHECKERBOARDING 

A common phenomenon in Galerkin finite element formulations known as "checker- 
boarding" (references 26, 27), or "cusping" (reference 21) is sometimes observed 
in areas where large mean flow gradients are present or where mesh resolution is 
inadequate. Sometimes the presence of a localized problem may cause an insta- 
bility which spreads and contaminates the entire solution. 

The presence of tolerable levels of contamination is shown in Figures 7.25 
through 7.27. These three figures are all parts of the same three AFE solu- 
tions discussed earlier. Figure 7.25 represents the axial variation of the 
axial components of acoustic velocity amplitude along a line close to the wall 
for case 3. Due to the presence of the boundary layer (where large mean vel- 
ocity gradients exist), small oscillations are superimposed on the true solu- 
tion. A magnified example from the same case for the radial velocity component 
is shown in Figure 7.26. Here the high level of contamination is tolerable due 
to the small overall magnitude of the radial velocity component compared to the 
axial component and its negligible contribution to the overall solution. 
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Even in apparently very stable solutions, some level of checkerboarding is fre- 
quently present. Figure 7.27 shows the effects of greatly magnifying the scale 
on the linear element solution radial plot (Figure 7.13) from the AFE no flow 
analysis of case 1. 
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Figure 7.5 Streamlines for Case 3 
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Figure 7.6 Mach Number Contours from .24 to .68 ( intervals 04) for Case 3 



Figure 7 7 Velocity Contours from 80 to 230 ( interval. 04) for Case 3 
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Figure 7.12 Axial Plot of Acoustic Flux Integrated over 
Successive Radial Planes for Case 1 (No Flow) 
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7.13 Radial Plot of Pressure Amplitude for Case 1 (No Flow), 
Axial Index 50, z = .3763m 


AXIAL INDEX SO Z = .3763 


-i 1 1 I i 




DUCT RADIUS 


Figure 7.14 Radial Plot of Pressure Phase for Case 1 (No Flow), 
Axial Index 50, z = .3763m 
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Figure 7.15 Axial Plot of the Real Part of Outer Wall Specific 
Acoustic Admittance for Case 1 (No Flow) 
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Figure 7.16 Axial Plot of the Imaginary Part of Outer Wall Specif 
Acoustic Admittance for Case 1 (No Flow) 
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8 SUMMARY 


This report describes an advanced system for the analysis of sound propagation 
in ax i symmetric (or planar) ducts containing near sonic compressible mean flows. 

The system derives its unique capabilities from the integration of a thoroughly 
tested computer code for derivation of the mean flow with an established finite 
element code for solution of the acoustic problem. 

Also of significance are steps which were taken to make the system "friendly" 
to the user. These include the addition of powerful interactive graphics for 
interpretation and analysis of results and the introduction of simplified 
operating procedures through use of advanced features of CDC-NOS on NASA 
Langley Research Center's Cyber computers. 

The Adam system is considered to represent a significant step forward in facil- 
itating advanced aeroacoustic duct analyses both for research and design. 
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APPENDIX 1 


Sheet 1 


Overall Input Data 


▼ 

input tape? 

output tape? 


2 

y 

T or F 

T or F 


1 STC 





Mach number, ambient pressure and temperature, fluid properties 
$A (1.) (1.) (1.) (1.4) 

MACHO , T0S= , PS0= , RG= , GAM= . 

Highlight radius, maximum body radius, body closure tolerance 

RHL= , RM= , TTE= , 

Axisymmetric or Planar? 

T or F 
AX I = , 

spatial grid refinement criteria, see notes 

GR ( 1 ) = , , , , , , , 

SGR(1)= , , , , , , , 

NGR= _, 

GZ(1)= , , , , , , , 

SGZ( 1 )= , , , , , , , 

NGZ= 


maximum Mach number increment between grid points 


streamwise 

direction 

(0.1) 


normal 

direction 

(0.1) 


VMG1= VMG2- 


maximum number curvature damping 

of refinements (1.) 


incompressible start 
( 0 ) 


MAXIT= RH0C= 


N0DENS= 


fluid reference temperature, reference viscosity Sutherland constant, unit 
conversion constant 

(518. 7°r) (10" 6 lbm/in.sec) (198.6° R) 

TREF= , MUREF= , SC0N= , 

$ 


74 



Sheet 2 



boundary 

name 

(1-6 characters) 

T ▼ T 


2 BDY 

upper boundary? angle input? 

T or F T-no, F-yes 



$ 


’ Coordinates 

channel 

name 

1-6 characters) 

T 

24 


boundary layer? equiv. flat plate dis- 

T-yes, (F-no) tance to boundary layer 

origin 

( 0 .) 

BL= , CAPX1= 
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Sheet 2a 


Pressure Boundaries (Optional) 

These data may be included as optional parameters between the $A and $ lines in 
Sheet 2. 

Total number of points on all pressure boundaries (MAX = 10) 

NZP = , 

Axial distance 

ZP(1) = , , , , 

Pressure or Velocities 

PPS(l) = , , , , 

(If two pressure boundaries exist, then data for the second follow data for 
first in above arrays). 


Number of above values applying to first boundary, P/V indicator ( = 0 for 
pressure, = 1 for velocity) 

NZP1 = , PSIV = (0) , 

Sheet 3 


Channel Flow Properties 

Channel name 

▼ ▼ ▼ 

2 4 14 

▼ ▼ ▼ 

3 CHN 


$A 


ratio of 
specific 
heats 
(1.4) 


gas 

constant 


(1.0) 


flow rate 
may be 
adjusted? 
(T) or F 


GAM- , RG= , VARY- 


stagnation properties 

total temp total pressure 

TT0- , PT0- , 

Mach no. static temp static pressure 

MACHO , TSO , PSO 

flow area normalized by A^ L 


AO- 

$ 
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Sheet 3a 

Restart Procedure for Merging STC 
Solution with Boundary Layer 

These data follow Sheet 3 and are simply repeats of Sheet 1 with 11 input tape" 
set to "true". To use this procedure, "output tape" must be set to "true" on 
the original Sheet 1. 

column 

▼ 

2 Input tape? Output tape 

T T or F T or F 

1 STC T T 

$A 

$ 

1 STC T T 

$A 

$ 

1 STC T T 

$A 

$ 

The above example performs three restarts. More may be added as required until 
it is seen that the solution is fully converged. 
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Specimen Card Input 


NAME= 

ADDRESS= 

IDENT=CONICAL-CIRCULAR CONSTRICTION A/A*=2.0351 

1 STC F T 

$A 

MACH0=0. 14,TS0=516.665,PS0=14. 5,RG=1716.0, 
RHL=5.895,RM=5.895,TTE=0.0, 

AXI=T, 

GR( 1 1=0.0,5.895, 

SGR ( 1 1=0.75,0.75, 

NGR=2, 

GZ(l)=-53. 0,-35. 0,-20. 0,0. 0,20. 0,35. 0,53.0, 

SGZ( 1) =3. 0,3. 0,0. 5, 0.1, 0.5, 3. 0,3.0, 

NGZ=7, 

VMG1=50.0,VMG2=75.0, 

MAXIT=2,PRPRN=-1.0, 

$ 

2 BDY CL INT 

$A 

UPPER=F,ZRONLY=F,BL=F, 

B(l)=-53. 0,0. 0,0.0, 

53.0,0.0,0.0, 

$ 

2 BDY CNLCIR INT 

$A 

UPPER=T,ZRONLY=F,BL=T, 

B( 1 ) =-53. 0,5. 895,0.0, 

-29.475, 5.895,0.0, 

-29.4687,5.8948, -3.6075, 

-3.0203, 4.2271, -3.6075, 

-2.9303, 4.2215, -3.5, 

-2.5121, 4.1978, -3.0, 

-2.0937, 4.1777 -2.5, 


-1.6752, 

4.1612, 

-2.0, 

-1.2565, 

4.1484, 

-1.5, 

-0.8377, 

4.1393, 

-1.0, 

-0.4189, 

4.1338, 

-0.5, 

0.0, 

4.132, 

0.0, 

0.4189, 

4.1338, 

0.5, 

0.8377, 

4.1393, 

1.0 

1.2565, 

4.1484, 

1.5 

1.6752, 

4.1612, 

2.0, 

2.0937, 

4.1777, 

2.5 

2.5121, 

4.1978, 

3.0 

2,9303, 

4.2215, 

3.5, 

3.0202, 

4.2271, 

3.6075, 

29,4687, 

5.8948, 

3.6075, 

29.475, 

5.895,0. 

. 0 , 

53.0,5.895,0.0, 



$ 
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3 CHN I NT 

$A 

RG= 1716.0, VARY=T , A0= 1.0, 

$ 

1 STC T T 

$A MAXIT=4, $ 

1 STC T T 

$A 

MACHO . 295 , TS0=509 . 8 1 5 , PS0= 1 3 . 8385 , 
MAXIT=7, 

1 STC T T 

$A MAXIT=7, $ 

1 STC T F 

$A MAXIT=8,PRPRN=1 .0,$ 
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Appendix 1 - Addendum - Additional STC P rogram Controls 


Part A - Optional input which may follow sheet 1 

A large number of input quantities that are additional to those described in the 
main body of this report may be specified. These items are not normally required 
for execution of the STC program and consist in general of controls for special 
program options and input data to modify preset or initialized constants or par- 
ameters. 

Only the following selected few of these optional input parameters are described 
here (for a full list see reference 17): 

Flow Balance and Flow Adjustment Tolerances 




Preset Value 

T0LES2 

Final solution tolerance on streamline 
position error (=MAXES2) 

.001 

TOLINR 

Intermediate solution tolerance on MAXES2 

.05 

NINNER (1) 

Number of inner iterations without grid 
refinement. Specify up to MAXIT Values 

16*10 

CNVF (1) 

Fractional percentage of total point 
movement to be used at a given refine- 
ment level. Specify up to MAXIT Values. 

10*1.0 

TOLWF 

Tolerance on fractional flow adjustment 
needed to meet trailing edge pressure 
closure condition. 

.001 


(A more detailed explanation of the function of these para- 
meters is contained in the following section). 

Grid Refinement Parameters 

The following items may be used to control the length of orthogonal 
lines and streamlines as the grid is refined. If set equal to zero, 
they ensure that all new grid lines will span the entire field. 

CRXSL Control for extending streamlines .375 

CRXOL Control for extending orthogonal lines .375 

across a subsonic region. 

Part B - Notes on interpretation of solution history output 

This output section describes the progress of the STC iterative solution 
algorithm. Here, the behavior of MAXES2 relative to LIMES2 indicates con- 
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vergence by a decreasing absolute value as the solution proceeds. Any erratic 
behavior in MAXES2 indicates convergence problems. 

Sometimes MAXES2 will start to converge and then show instability. By identi- 
fying the appropriate NINNER, the solution can be forced to the previous level 
of grid refinement by rerunning the problem with the NINNER control inserted. 
Convergence difficulties for this type of problem may sometimes be overcome by 
establishing a refined grid at a lower Mach number and then slowly raising the 
Mach number to the desired level using restart cases. 

Definitions of parameters listed in the solution history follow: 


NREFINE - The number of times the grid has been refined. The maximum number of 
grid refinements is controlled by input parameter MAX IT. 

GRIDPTS - The number of points in the flow field. 

INRCTR - The number of inner iterations at the current grid refinement level. 

These are controlled by built in tolerances on the flow balance, but 
may be overidden by inputting NINNER(M) where M is the level of grid 
refinement. 


NSSPTS - Number of imbedded supersonic points in the flow field. 

NSWEEPS - Number of sweeps through the matrix solution on the current inner 
iteration. 


MAXDS2 - Maximum streamline adjustment actually made during the current inner 
iteration. 

LIMES2 - TOLINR * (current grid size), where TOLINR is an optional input para- 
meter specifying a tolerance on grid development, i.e., grid devel- 
opment in the inner loop continues provided MAXES2 < LIMES2 (unless 
preceeded by other conditions). 

(r,z) - Field coordinates of point with maximum streamline movement. 


The utility of TOLINR is that it permits a looser tolerance on inner iterations 
while grid development is taking place. Finally, when grid development is 
satisfactorily achieved, then a final tight tolerance (T0LES2) assumes control, 
and flow field definition continues until MAXES2 < T0LES2 * CLEN where CLEN = 
(Average of SGR's and SGZ's). 
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APPENDIX 2 


Format of AFE12 Output File 

(Total number of Records = 2N + 4) 

Record 1 - AFE Program Input Data Identifiers 

M - Number of elements in radial direction 
N = Number of elements in axial direction 
NPHPAR = Number of acoustic physical parameters at a node 
NNDERIV = Number of nodal derivatives of each parameter 
FLOW = "True" if flow data on GE file was used 
= "False" for no flow 
FREQ = Acoustic Frequency 


Record 2 - STC Program Input Data Identifiers 

NAME (6), IDENT (6), MACHO, TSO, PSO, RG, 
NAMEBDY (2), UPPER (2), BL (2) 


Record 3 - Aerodynamic parameters from STC program interpolated on AFE mesh 
AERODAT (Ml, Nl, 2) = Mach number and flow angle at each node 

Record (3+J) - Data for AFE row #J (i.e. nodal rows J and J+l) 

(N Records of length 4 (M+l) + 16M) 

RZ (2, Ml, 2) = Nodal coordinates 

Index 3=1 : J'th nodal row 

= 2 : (J+l ) ' t h nodal row 

Index 2 = 1 to Ml : Radial node number 

Index 1=1 : R coordinates 

= 2 : Z coordinates 

JAC (4, 4, M) = Jacobian matrices for J'th AFE row 

Index 3 = 1 to M for elements in radial direction 

Indices 1&2 : e.g. = J^ = JAC (1, k, I) 

J 4k = JAC (4, k, I) 
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Record (3 + N + J) - Solution vectors for J'th radial 

(N+l) Records of length NPHPAR * NNDERIV * (M+l)* 2+1 
parts of solution vectors with flux in last location. 

Order of variables is as follows: 

9U 9u 9 2 U 9 V 9 V 9 2 V 9W 9W 9 2 W 

* 3r’ 9z ’ 9r9z’ ’ ar ’ az ’ araz ’ ’ ar ’ az ’ ar az ’ P’ 


nodal rows (and flux). 
= real and imaginary 

9p 9p 9 2 p 
ar’ az ’ 9r 9z * 


83 


r 


1. Report No. 2. Government Accession No. 

NASA CR-3668 


4. Title and Subtitle 

ADAM - AN AXISYMMETRIC DUCT AEROACOUSTIC 
MODELING SYSTEM 


7. Author(s| 


3. Recipient’s Catalog No. 


5. Report Date 

January 1983 


6. Performing Organization Code 


A. Louis Abrahamson 


8. Performing Organization Report No. 


10. Work Unit No. 


11. Contract or Grant No. 


13. Type of Report and Period Covered 

Contractor Report 


14. Sponsoring Agency Code 

505-31 33 1 3 
L-34073B 


9. Performing Organization Name and Address 

COMTEK 

123 East Woodland Road 
Grafton, Virginia 23692 


12. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 
Washington, D.C. 20546 

15. Supplementary Notes 

Langley Technical Monitor: Harold C. Lester 
Final Report 

16. Abstract 

ADAM is an acronym for Axisymmetric Duct Aeroacoustic Modeling and is the name of an 
interconnected system of computer programs for analysing the propagation and attenuation 
of sound in aeroengine ducts containing realistic compressible subsonic mean flows. It 
has been developed primarily for research directed towards the reduction of noise emitted 
from turbofan aircraft engines and includes as two basic components a streamtube curvature 
program for determination of the mean flow and a finite element code for solution of the 
acoustic propagation problem. 

The system, which has been specifically tailored for ease of use, is presently installed 
at NASA Langley Research Center on a Control Data Cyber 175 Computer under the 
NOS Operating System employing a Tektronix terminal for interactive graphics. 

This report describes the scope and organization of the ADAM System and includes a 
users guide, examples of input data and results for selected cases. 


17. Key Words (Suggested by Author(s)) 

Duct Acoustics 
Sound Propagation 
Noise Control 


18. Distribution Statement 


19. Security Classif. (of this report! 

Unclassified 


Unclassified - Unlimited 


Subject Category 71 


20. Security Classif. (of this page) 

Unclassified 


21. No. of Pages 22. Price 



For sale by the National Technical Information Service, Springfield, Virginia 22161 


NASA-Langley, 1983 







